#include <stdio.h>

#define TAMANHO_BLOCO 10

/* Estutura de dados de uma matriz quadrada (MatQ) */
struct MatQ{
  int ordem; /* ordem da matriz quadrada */
  float *elementos; /* elementos da matriz (salvos na forma de um vetor) */
  int tamanhoLinha; /* tamanho da linha */
};
typedef struct MatQ MatQ;

/* Estutura de dados de um vetor */
struct Vet{
  int tamanho; /* tamanho do vetor */
  float *elementos; /* elementos do vetor */
};
typedef struct Vet Vet;

/*
 * Cria uma matriz quadrada (MatQ) na memoria RAM (no PC)
 * parametros:
 *  mat: matriz a ser criada
 *  ordem: ordem da matriz
 */
void criarMatQPC(MatQ *mat, int ordem){
  mat->ordem = mat->tamanhoLinha = ordem; /* define a ordem da matriz */
  
  mat->elementos = (float *)malloc(sizeof(float) * ordem * ordem); /* aloca a memoria para os elementos */
}

/*
 * Cria uma matriz quadrada (MatQ) na memoria da GPU
 * parametros:
 *  mat: matriz a ser criada
 *  ordem: ordem da matriz
 */
void criarMatQGPU(MatQ *mat, int ordem){
  mat->ordem = mat->tamanhoLinha = ordem; /* define a ordem da matriz */
  
  cudaMalloc(&(mat->elementos), sizeof(float) * ordem * ordem); /* aloca a memoria para os elementos */
}

/*
 * Copia uma matriz quadrada (MatQ) da memoria da GPU para memoria RAM (matPC = matGPU)
 * parametros:
 *  matPC: matriz na memoria ram, que recebera a copia
 *  matGPU: matriz na GPU, que sera copiada
 */
void copiarMatQParaPC(MatQ *matPC, MatQ *matGPU){
  /* faz a copia da GPU para a RAM */
  cudaMemcpy(matPC->elementos, 
              matGPU->elementos, 
              sizeof(float) * matPC->ordem * matPC->ordem, 
              cudaMemcpyDeviceToHost);
}

/*
 * Copia uma matriz quadrada (MatQ) da memoria RAM para memoria da GPU (matGPU = matPC)
 * parametros:
 *  matGPU: matriz na GPU, que recebera a copia
 *  matPC: matriz na memoria ram, que sera copiada
 */
void copiarMatQParaGPU(MatQ *matGPU, MatQ *matPC){
  /* faz a copia da RAM para a GPU */
  cudaMemcpy(matGPU->elementos, 
              matPC->elementos, 
              sizeof(float) * matGPU->ordem * matGPU->ordem, 
              cudaMemcpyHostToDevice);
}

/* 
 * Desaloca uma matriz da memoria RAM 
 * parametros:
 *  mat: matriz a ser desalocada
 */
void deletarMatQPC(MatQ *mat){
  free(mat->elementos);
}

/* 
 * Desaloca uma matriz da memoria da GPU
 * parametros:
 *  mat: matriz a ser desalocada
 */
void deletarMatQGPU(MatQ *mat){
  cudaFree(mat->elementos);
}

/*
 * Cria um vetor na memoria RAM (no PC)
 * parametros:
 *  v: vetor a ser criado
 *  tamanho: tamanho do vetor
 */
void criarVetPC(Vet *v, int tamanho){
  v->tamanho = tamanho; /* define o tamanho */
  
  v->elementos = (float *)malloc(sizeof(float) * tamanho); /* aloca os elementos */
}

/*
 * Cria um vetor na memoria da GPU
 * parametros:
 *  v: vetor a ser criado
 *  tamanho: tamanho do vetor
 */
void criarVetGPU(Vet *v, int tamanho){
  v->tamanho = tamanho; /* define o tamanho */
  
  cudaMalloc(&(v->elementos), sizeof(float) * tamanho); /* aloca os elementos */
}

/*
 * Copia um vetor da memoria da GPU para memoria RAM (vPC = vGPU)
 * parametros:
 *  vPC: vetor na memoria ram, que recebera a copia
 *  vGPU: vetor na GPU, que sera copiado
 */
void copiarVetParaPC(Vet *vPC, Vet *vGPU){
  /* faz a copia da GPU para a RAM */
  cudaMemcpy(vPC->elementos, 
              vGPU->elementos, 
              sizeof(float) * vPC->tamanho, 
              cudaMemcpyDeviceToHost);
}

/*
 * Copia um vetor da memoria RAM para a memoria da GPU (vGPU = vPC)
 * parametros:
 *  vGPU: vetor na GPU, que recebera a copia
 *  vPC: vetor na memoria ram, que sera copiado
 */
void copiarVetParaGPU(Vet *vGPU, Vet *vPC){
  /* faz a copia da RAM para a GPU */
  cudaMemcpy(vGPU->elementos, 
              vPC->elementos, 
              sizeof(float) * vGPU->tamanho, 
              cudaMemcpyHostToDevice);
}

/*
 * Copia v2 para v1 (os dois vetores estao na GPU)
 * parametros:
 *  v1: vetor na GPU, que recebera a copia
 *  v2: vetor na GPU, que sera copiado
 */
void copiarVetGPU(Vet *v1, Vet *v2){
  /* copia os elementos de v1 para v2*/
  cudaMemcpy(v1->elementos, 
              v2->elementos, 
              sizeof(float) * v2->tamanho, 
              cudaMemcpyDeviceToDevice);
}

/* 
 * Desaloca um vetor da memoria RAM
 * parametros:
 *  v: vetor a ser desalocado
 */
void deletarVetPC(Vet *v){
  free(v->elementos);
}

/* 
 * Desaloca um vetor da memoria da GPU
 * parametros:
 *  v: vetor a ser desalocado
 */
void deletarVetGPU(Vet *v){
  cudaFree(v->elementos);
}

/* 
 * Troca os elementos de dois vetores
 * parametros:
 *  vet1: ficara com os elementos de vet2
 *  vet2: ficara com os elementos de vet1
 */
void trocarVetsGPU(Vet *vet1, Vet *vet2){
  float *aux = vet1->elementos; /* ponteiro auxiliar */
  
  /* troca os ponteiros dos elementos dos vetores */
  vet1->elementos = vet2->elementos;
  vet2->elementos = aux;
}

/* 
 * Retorna um elemento de uma matriz que esta na memoria RAM
 * parametros:
 *  mat: matriz onde esta o elemento
 *  linha: linha do elemento
 *  coluna: coluna do elemento
 * retorno:
 *  elemento da matriz
 */
float getElemento(const MatQ *mat, int linha, int coluna){
  return mat->elementos[linha * mat->tamanhoLinha + coluna]; /* posicao do elemento no vetor */
}

/* 
 * Define o valor de um elemento de uma matriz que esta na memoria RAM
 * parametros:
 *  mat: matriz onde esta o elemento
 *  linha: linha do elemento
 *  coluna: coluna do elemento
 *  valor: valor do elemento
 */
void setElemento(MatQ *mat, int linha, int coluna, float valor) {
  mat->elementos[linha * mat->tamanhoLinha + coluna] = valor; /* posicao do elemento no vetor */
}

/* 
 * Retorna um elemento de uma matriz que esta na GPU
 * parametros:
 *  mat: matriz onde esta o elemento
 *  linha: linha do elemento
 *  coluna: coluna do elemento
 * retorno:
 *  elemento da matriz
 */
__device__ float getElementoGPU(const MatQ *mat, int linha, int coluna){
  return mat->elementos[linha * mat->tamanhoLinha + coluna]; /* posicao do elemento no vetor */
}

/* 
 * Define o valor de um elemento de uma matriz que esta na GPU
 * parametros:
 *  mat: matriz onde esta o elemento
 *  linha: linha do elemento
 *  coluna: coluna do elemento
 *  valor: valor do elemento
 */
__device__ void setElementoGPU(MatQ *mat, int linha, int coluna, float valor) {
  mat->elementos[linha * mat->tamanhoLinha + coluna] = valor; /* posicao do elemento no vetor */
}

/* 
 * Imprime uma matriz que esta na RAM
 * parametros:
 *  mat: matriz a ser impressa
 */
void imprimir(MatQ *mat){
  int i, j, ordem = mat->ordem;
  
  for(i = 0; i < ordem; ++i){
    for(j = 0; j < ordem; ++j){
      printf("%f ", getElemento(mat, i, j));
    }
    printf("\n");
  }
}

/* 
 * Calcula uma iteracao da aproximacao da solucao do sistema
 * parametros:
 *  matA: matriz A (coeficientes das variaveis nas equacoes do sistema)
 *  vetB: vetor B
 *  x: valor da aproximacao do vetor X nessa iteracao (retorno)
 *  xAnt: valor da aproximacao do vetor X na iteracao anterior
 */
__global__ void calcularIteracao(MatQ matA, Vet vetB, Vet x, Vet xAnt) {
  int linhaBloco = blockIdx.y; /* linha do bloco de threads */
  
  int linhaThread = threadIdx.y; /* linha da thread no bloco */
  int colunaThread = threadIdx.x; /* coluna da thread no bloco */
  
  int linhaMat = linhaBloco * TAMANHO_BLOCO + linhaThread; /* linha da matriz A que a thread atual ira acessar */
  int colunaMat = colunaThread; /* coluna da matriz A que a thread atual ira acessar */
  
  float aux = 0.0; /* variavel auxiliar */

  /* se as variaveis linhaMat e colunaMat estiverem nos limites da matriz */
  if(linhaMat < matA.ordem && colunaMat < matA.ordem) {
    for (int m = 0; m < (matA.ordem / TAMANHO_BLOCO); ++m) {
      __shared__ float cacheA[TAMANHO_BLOCO][TAMANHO_BLOCO]; /* chache de parte da matriz A para a memoria compartilhada das threads */
      __shared__ float cacheB[TAMANHO_BLOCO]; /* chache de parte do vetor B para a memoria compartilhada das threads */
      
      cacheA[linhaThread][colunaThread] = getElementoGPU(&matA, linhaMat, colunaMat); /* copia o elemento da matriz para o cache */
      
      /* se a thread for da linha zero copia o elemento do vetor B (para evitar acessos de gravacao concorrentes ao cache do vetor B) */
      if(linhaThread == 0){
        cacheB[colunaThread] = xAnt.elementos[colunaMat]; /* copia o elemento para o cache */
      }
      
      __syncthreads(); /* sincroniza as threads */
      
      /* multiplica cada elemento do cache da matriz pelo seu correspondente no vetor */
      cacheA[linhaThread][colunaThread] *= cacheB[colunaThread];
      
      __syncthreads(); /*sincroniza as threads */
      
      /* se a thread for a primeira da linha ela realizara a soma de A[i][j]X[j] */
      if(colunaThread == 0){
        for (int e = 0; e < TAMANHO_BLOCO; ++e){
          aux += cacheA[linhaThread][e];
        }
      }
        
      colunaMat += TAMANHO_BLOCO; /* incrementa o contador de colunas */
      
      __syncthreads(); /* sincroniza as threads */
    }

    /* se a thread for a primeira da linha ela salva o x no vetor X (B[i] - Soma(A[i][j]X[j])) */
    if(colunaThread == 0){
      x.elementos[linhaMat] = vetB.elementos[linhaMat] - aux;
    }
    
    __syncthreads(); /* sincroniza as threads */
  }
}

/* 
 * Normaliza a matriz A e o vetor B
 * parametros:
 *  matA: matriz A
 *  matANorm: matriz A normalizada (retorno)
 *  vetB: vetor B
 *  vetBNorm: vetor B normalizado (retorno)
 */
__global__ void normalizar(MatQ matA, MatQ matANorm, Vet vetB, Vet vetBNorm){
  int linhaThread = threadIdx.y; /* linha da thread no bloco */
  int colunaThread = threadIdx.x; /* coluna da thread no bloco */
  int linhaBlocoThread = blockIdx.y; /* linha do bloco de threads */
  int colunaBlocoThread = blockIdx.x; /* coluna do bloco de threads */

  int linhaMat = linhaBlocoThread * TAMANHO_BLOCO + linhaThread; /* linha da matriz A correspondente a thread */
  int colunaMat = colunaBlocoThread * TAMANHO_BLOCO + colunaThread; /* coluna da matriz A correspondente a thread */
  
  __shared__ float diagonal[TAMANHO_BLOCO]; /* cache da diagonal da matriz A */
  
  /* se as variaveis linhaMat e colunaMat estiverem nos limites da matriz */
  if(linhaMat < matA.ordem && colunaMat < matA.ordem) {
    /* se linha e a coluna da thread sao iguais ela deve salvar no cache um elemento da diagonal */
    if(linhaThread == colunaThread){
        diagonal[linhaThread] = getElementoGPU(&matA, linhaMat , linhaMat); /* salva o elemento correspondente a sua linha */
    }
    
    __syncthreads(); /* sincroniza as threads */
    
    /* se a thread corresponde a um elemento da diagonal da matriz */
    if(linhaMat == colunaMat){
      vetBNorm.elementos[linhaMat] = vetB.elementos[linhaMat] / diagonal[linhaThread]; /* normaliza um elemento de B */
      setElementoGPU(&matANorm, linhaMat, colunaMat, 0); /* zera um elemento da diagonal da matriz A */
    }else{
      float aux = getElementoGPU(&matA, linhaMat, colunaMat);
      setElementoGPU(&matANorm, linhaMat, colunaMat, aux / diagonal[linhaThread]); /* normaliza um elemento */
    }
  }
}

/* 
 * Calcula o erro atual
 * parametros:
 *  x: x na iteracao atual
 *  xAntes: x na iteracao passada
 * retorno:
 *  erro na iteracao atual
 */
float calcularErro(Vet x, Vet xAntes){
  float absDifAux, /* variavel auxiliar para o calculo da maior diferenca absoluta */
        absDifMaior = fabs(x.elementos[0] - xAntes.elementos[0]), /* maior diferenca absoluta */ 
        absXMaior = fabs(x.elementos[0]), /* maior x absoluto */
        absXAux; /* variavel auxiliar paar o calculo da maior x absoluto */
  int i, /* contador */ 
      n = x.tamanho; /* tamanho dos vetores */
  
  for(i = 1; i < n; ++i){
    absXAux = fabs(x.elementos[i]); /* calcula o x absoluto atual */
    absDifAux = fabs(x.elementos[i] - xAntes.elementos[i]); /* calcula a diferenca absoluta atual */
    
    /* verifica se o x absoluto atual e o maior ate o momento */
    if(absXAux > absXMaior){
      absXMaior = absXAux; /* armazena o maior x absoluto */
    }
    
    /* verifica se a diferenca absoluta atual e a maior ate o momento */
    if(absDifAux > absDifMaior){
      absDifMaior = absDifAux; /* armazena a maior diferenca absoluta */
    }
  }
  
  /* retorna o erro */
  return (absDifMaior / absXMaior);
}

/* 
 * Calcula a aproximacao da solucao do sistema
 * parametros:
 *  matA: matriz A
 *  vetB: vetor B
 *  nIteracoes: numero maximo de iteracoes a serem executadas para o calculo da aproximacao (condicao de parada)
 *  linhaAvaliada: linha da matriz A (equacao do sistema) a ser usada para avaliar a aproximacao
 *  erroMaximo: maior erro permitido (condicao de parada)
 */
void resolverSistema(MatQ *matA, Vet *vetB, int nIteracoes, int linhaAvaliada, float erroMaximo){
  int i, /* contador */
      ordemMat = matA->ordem; /* ordem da matriz A */
  MatQ matAGPU, /* copia da matriz A na GPU */
       matANormGPU; /* matriz A normalizada na GPU (A*) */
  Vet vetBGPU, /* copia do vetor B na GPU */
      vetBNormGPU, /* vetor B normalizado na GPU (B*) */
      vetXGPU, /* aproximacao de X na iteracao atual armazenado na GPU */
      vetXAntGPU, /* aproximacao de X na iteracao anterior armazenado na GPU */
      resultados, /* vetor resultado da aproximacao (ao fim de todas a itercoes) armazenado na RAM */
      vetX, /* aproximacao de X na iteracao atual armazenado na RAM */
      vetXAnt; /* aproximacao de X na iteracao anterior armazenado na RAM */
  float erro; /* erro na iteracao atual */
  
  criarMatQGPU(&matAGPU, ordemMat); /* aloca a matriz A na GPU */
  copiarMatQParaGPU(&matAGPU, matA); /* copia a matriz A da RAM para a GPU */
  
  criarMatQGPU(&matANormGPU, ordemMat); /* aloca a matriz A normalizada (A*) na GPU */
  
  criarVetGPU(&vetBGPU, ordemMat); /* aloca o vetor B na GPU */
  copiarVetParaGPU(&vetBGPU, vetB); /* copia o vetor B da RAM para a GPU */
  
  criarVetGPU(&vetBNormGPU, ordemMat); /* aloca o vetor B normalizado (B*) na GPU */
  
  /* define as dimensoes dos blocos e do grid usados na normalizacao */
  dim3 dimBlock(TAMANHO_BLOCO, TAMANHO_BLOCO);
  dim3 dimGrid(ordemMat / dimBlock.x + 1, ordemMat / dimBlock.y + 1);
  
  /* normaliza a matriz A e o vetor B (matAGPU e vetBGPU) e armazena os resuldados nas variaves matANormGPU e vetBNormGPU */
  normalizar<<<dimGrid, dimBlock>>>(matAGPU, matANormGPU, vetBGPU, vetBNormGPU);
  
  criarVetGPU(&vetXGPU, ordemMat); /* aloca o vetor X na GPU */
  copiarVetGPU(&vetXGPU, &vetBGPU); /* copia o vetor X da RAM para a GPU */
  
  criarVetGPU(&vetXAntGPU, ordemMat); /* aloca o vetor X da iteracao anterior na GPU */
  copiarVetGPU(&vetXAntGPU, &vetBNormGPU); /* copia o vetor X da iteracao anterior da RAM para a GPU */
 
  /* define as dimensoes do grid usado para resolver as iteracoes do sistema */
  dim3 dimGrid2(1, ordemMat / dimBlock.y);

  criarVetPC(&resultados, ordemMat); /* aloca o vetor resultados na RAM */
  
  criarVetPC(&vetX, ordemMat); /* aloca o vetor X na RAM */
  criarVetPC(&vetXAnt, ordemMat); /* aloca o vetor X da iteracao anterior na GPU */
  
  erro = 10; /* inicia o erro com um valor grande para que o loop entre na primeira iteracao */
  /* calcula as iteracoes ate que o erro esteja em um valor permitido ou que o numero maximo de iteracoes seja alcancado */
  for(i = 0; i < nIteracoes && erro > erroMaximo; ++i){
      
    /* calcula o vetor X para a iteracao atual */
    calcularIteracao<<<dimGrid2, dimBlock>>>(matANormGPU, vetBNormGPU, vetXGPU, vetXAntGPU); 
    
    /* copia os vetores X das ultimas iteracoes para a memoria */
    copiarVetParaPC(&vetX, &vetXGPU);
    copiarVetParaPC(&vetXAnt, &vetXAntGPU);
    
    /* calcula o erro na iteracao atual */
    erro = calcularErro(vetX, vetXAnt);

    /* troca os vetores dos x e x da iteracao anterior para que possa executar a proxima iteracao */
    trocarVetsGPU(&vetXAntGPU, &vetXGPU);
  }
  
  /* copia a ultima aproximacao para o vetor resultados */
  copiarVetParaPC(&resultados, &vetXAntGPU);

  printf("Iterations: %d\n", i); /* imprime o numero de iteracoes */
  
  /* calcula o valor da linha a ser avaliada */
  float aux = 0; /* zera a variavel auxiliar que ira computar a soma */
  for(i = 0; i < ordemMat; ++i){ /* computa a soma da linha */
    aux += resultados.elementos[i] * getElemento(matA, linhaAvaliada, i);
  }
  printf("RowTest: %d => [%f] =? %f\n", linhaAvaliada, aux, vetB->elementos[linhaAvaliada]);  /* imprime o resultado da soma da linha e o valor esperado */
    
  /* desaloca as variveis alocadas na RAM e na GPU */
  deletarMatQGPU(&matAGPU);
  deletarVetGPU(&vetBGPU);
  
  deletarMatQGPU(&matANormGPU);
  deletarVetGPU(&vetBNormGPU);
  
  deletarVetGPU(&vetXGPU);
  deletarVetGPU(&vetXAntGPU);
  deletarVetGPU(&vetXAntGPU);
  
  deletarVetPC(&vetX);
  deletarVetPC(&vetXAnt);
  deletarVetPC(&resultados);
}

int main(void){
  MatQ matA; /* matriz A */
  Vet vetB; /* vetor B */
  int ordemMat, /* ordem da matriz A */
      filaAvaliada, /* linha do sistema a ser usada na avaliacao */
      maxIteracoes, /* maximo de iteracoes para o calculo da aproximacao */
      i, /* contador */
      j; /* contador */
  float erroPermitido, /* erro maximo permitido */
        aux; /* variavel auxiliar usada na leitura da matriz A */
  
  /* le os parametros do sistema */
  scanf("%d %d %f %d", &ordemMat, &filaAvaliada, &erroPermitido, &maxIteracoes);

  /* aloca a matriz A na RAM */
  criarMatQPC(&matA, ordemMat);
  
  /* le os elementos da matriz A */
  for(i = 0; i < ordemMat; ++i){ 
    for(j = 0; j < ordemMat; ++j){
      scanf("%f", &aux);
      setElemento(&matA, i, j, aux);
    }
  }

  /* aloca o vetor B na RAM */
  criarVetPC(&vetB, ordemMat);
  
  /* le os elementos do vetor B */
  for(i = 0; i < ordemMat; ++i){
    scanf("%f", &vetB.elementos[i]);
  }

  /* resolve o sistema (calcula a aproximacao da solucao) */
  resolverSistema(&matA, &vetB, maxIteracoes, filaAvaliada, erroPermitido);
  
  /* deleta a matriz A da RAM */
  deletarMatQPC(&matA);
  
  /* deleta o vetor B da RAM */
  deletarVetPC(&vetB);

  return 0;
}

